EPJ manuscript No. 

(will be inserted by the editor) 



Ab Initio Lattice Dynamics and Elastic Constants of ZrC 
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Abstract. Ab initio calculations and a direct method are applied to derive the phonon dispersion relations 
and phonon density of states for the ZrC crystal. The results are in good agreement with neutron scat- 
tering data. The force constants are determined from the Hellmann-Feynman forces induced by atomic 
displacements in a 2x2x2 supercell. The elastic constants are found using the deformation method and 
successfully compare with experimental data. 

PACS. 63.20.-e Phonons in crystal lattices - 71. 15. Mb Density functional theory 



Studies of transition-metal carbides have been carried 
out several times. Phenomenological model calculations [TJ 
managed to reproduce the phonon dispersion relations 
quite well, but the interactions were quite sophisticated , 
some of them including three-body interactions and dou- 
ble shell models. Recently, there have been several studies 
based on ab initio methods 4,5,6,7 . The electronic struc- 
ture, bulk modulus, and elastic constants, as well as the 
phonon dispersion relations have been found for TiC, TiN 
and TiO compounds by means of first-principles total- 
energy calculations JHIE] - Generally, the calculated values 
show good agreement with experiments. This provides a 
motivation to extend the investigation to heavier transi- 
tion metals and to treat zirconium carbide, ZrC. This crys- 
tal is an important material in nuclear energy technology, 
since it provides a filling medium for fuel particles. It is 
also used for surface hardening and coverage of cathodes 
of X-ray sources. In this paper we continue to study the 
lattice dynamics and elastic properties of transition-metal 
compounds. See [5] and references given there for a more 
detailed description of the method. 

The transition-metal carbide compounds, of which ZrC 
is a representative, are of considerable scientific and tech- 
nological interest because of their striking mechanical prop- 
erties, extreme hardness combined with metallic electri- 
cal and thermal conductivities. The ZrC crystal has NaCl 
structure and it is usually non-stoichiometric, mainly ow- 
ing to carbon- vacancy defects. The phonon dispersion re- 
lations of ZrC have been measured along main symme- 
try directions by Smith et al pUlfTT] . Elastic constants, in 
turn, have been measured for two concentrations of carbon 
(ZrCo.89 and ZrCo.94) by Chang and Graham |12) . 

In this paper we extend the first-principle calculations 
to describe the phonon dispersion curves, phonon density 
and elastic constants of ZrC. The method which is used, 
is based on the total energy calculation and Hellmann- 
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Feynman (HF) forces. The phonon dispersion relations 
are calculated by the direct method and the Phonon pro- 
gram |131I14II15HT^1[T71[T5] . in which the force constants of 
the dynamical matrix are derived from HF forces. Alter- 
natively, one could use the linear-response method |19[ 
120] for the evaluation of ab initio phonon frequencies at 
a predetermined set of Brillouin zone points. Elastic con- 
stants and bulk modulus arc estimated by straightforward 
evaluation of energy derivatives with respect to the defor- 
mations. 

The energies and HF forces of ZrC crystal are calcu- 
lated by the method of total energy minimization, using 
norm-conserving pseudopotentials as an approximation for 
the atomic core-valence electron interaction [2T1I221I23] . 
This method allows to include the pressure in calculations 
as well as anharmonic effects. For the lattice dynamics cal- 
culations a 2x2x2 supercell with periodic boundary con- 
ditions and 64 atoms was used. For optimizing the struc- 
ture and for the direct calculation of stress-strain relations 
a 1 x 1 x 1 supercell was utilized. The ab initio total energy 
calculations were done with the CASTEP package [23] 
and standard pseudopotentials constructed within LDA 
approximation and provided within this package. The Zr 
pseudopotential treats 4d electrons as belonging to the va- 
lence band. All pseudopotentials were parametrized in the 
reciprocal space with 680 eV cut-off energy. 

Tests, which were made with the lxlxl supercell and 
with cut-off energies of 680 eV and 900 eV, showed that 
in the first case the cohesive energy is only 0.01 eV higher, 
and the lattice constant changes by only 0.0007 A(0.01%). 
Therefore, we used the 680 eV cut-off energy for all re- 
maining calculations. 

The Local Density Approximation (LDA) as well as the 
Generalized Gradient Approximation (GGA) were used 
for the exchange energy term of the valence states of the 
Hamiltonian [25j . The integration over the Brillouin zone 
was performed with weighted summation over wave vec- 
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tors generated by the Monkhorst-Pack scheme using 
grid spacings from 0.1 A -1 to 0.04 A -1 which lead to 
sets of k- vectors containing from 4 to 118 wave vectors 
for the lxlxl supercell and from 1 to 14 wave vectors 
for the 2x2x2 supercell. The convergence of the force 
constants was achieved for grid spacing below 0.05 A" 1 . 
Hence, we carried on the optimization for a 0.05 A" 1 grid 
spacing, which leads to 32 and 4 k-points for the lxlxl 
and the 2x2x2 supercells, respectively. However, we used 
a 64 wave- vector set for the LDA calculation of the elastic 
constants. The metallic character of the TiC compound, 
studied in [5] , indicated that the use of smearing does not 
substantially improve the quality of phonon dispersion re- 
lations. Thus, we did not use smearing in the case of ZrC. 
The use of non-metallic summation of the Brillouin zone 
is further justified by the fact that the electron density 
of states at the Fermi level is small for ZrC and similar 
compounds: TiC and HfC [Tll2"T]. and the related errors 
are of the same order as the errors of the DFT approach 
itself [8 . Namely, they are less then 0.05% for the lattice 
constant, less then 0.3% for the bulk modulus, and less 
then 5% for phonon frequencies. 
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Fig. 2. Phonon density of states g{uf) of the ZrC crystal and 
partial phonon density of states gi ta (u), where a — C or Zr 
and i — x, y or z are the displacement directions. These are 
the GGA results. Experimental data taken from Ref. 1 28 . 
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Fig. 1. Phonon dispersion relations of ZrC crystal calculated 
with a 2x2x2 supercell within GGA (full line) and LDA 
(dashed line) approximations. The experimental points are 
taken from Ref. |10llllj . 



The minimization of the total energy leads to the equi- 
librium lattice constant of the stoichiometric ZrC a = 
4.691A for LDA and a = 4.695A for GGA. These values 
could be compared with the experimental values of a = 
4.6994A for ZrC . 9 4 and a = 4.7004A for ZrCo.sgDS- The 
phonon dispersion relations, which correspond to temper- 
ature T=0, are shown in Fig. [T] These are compared with 
experimental phonon frequencies measured by inelastic 
neutron scattering [ilJlITT] . For GGA, the calculated force 
constants for all 10 coordination spheres of the 2x2x2 
supercell are presented in Table [1] 



From Fig.[T]one sees that the experimental data at the 
r point are about Auj — 0.05a; lower than the calculated 
values. A reason for this could be that the experimen- 
tal values are measured at room temperature and that 
anharmonic effects diminish the phonon frequencies. An- 
other reason may relate to the non-stoichiometric concen- 
tration of carbon, which diminishes the average phonon 
frequencies as well. In other ab-initio calculations for al- 
kali metals [T3], graphite and diamond [13] the experimen- 
tal phonon points also are at frequencies lower than cal- 
culated ones. On the same Fig. [T] we see the influence of 
the different exchange-correlation energy approximations 
on phonon frequencies. In this case, the LDA appears to 
be far less accurate than the GGA. 

By sampling the dynamical matrix at many wave vec- 
tors, one can calculate the phonon density of states g(uj), 
and the partial phonon density of states g x ,Zr(^>) and 
g Xt c{oS). The g(uj) describes the number of phonon fre- 
quencies in an interval around w, while g x ,Zr and g x ^c 
specify the number of phonon frequencies in an inter- 
val around w, but only taking into account vibrations of 
Zr and C atoms, respectively. The density-of-state func- 
tions calculated from the GGA data are shown in Fig. [5] 
They are conventionally normalized to j g(u>)duj — 1 and 
J g Xt01 (u))du> = 1/6, where a = Zr or C. The curves show 
that motions within acoustic dispersion curves are almost 
entirely due to Zr atoms. The sum of the zirconium density 
of states (g x ,Zr(^) + gy ,Zr(^) + gz,Zr(oj)) fits the total den- 
sity of states g(u>) below 10 THz. The part of g(ui) above 
the gap is mainly due to vibrations of C atoms. Thus, ZrC 
forms quite a special crystal, in which the heavy Zr atoms, 
form a frame for the elastic motion, and the C atoms vi- 
brate within the optical modes. We see in Table [TJ that 
the magnitudes of the force constants between Zr-C and 
between Zr-Zr and C-C at a similar distance remain of 
the same order. Thus, ZrC does not consist of two weakly 
bounded subsystems. The situation is similar in TiC [S]. 
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Table 1. Values of the non-zero elements of the force-constant matrices derived from a 2x2x2 supercell in GGA (in N/m). All 
force-constant matrices are symmetric. 



In Fig. [5] the experimentally determined phonon density 
of states taken from Ref. [35] is also shown. The agreement 
within the acoustic region is quite good. In the optical re- 
gion the experimental resolution of ~ 2.5 THz and the 
non-stoichiometry of carbon lead to considerable broad- 
ening. 
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Table 2. Experimental and calculated bulk modulus B and 
elastic constants Cij of the ZrC crystal, in units of 10 11 Nm -2 = 
IMBar. Values in parentheses are calculated from B = |(cii + 
2ci 2 ). 



The bulk modulus was calculated in the lxlxl su- 
percell from the relation of pressure to the lattice con- 
stant. These data were fitted to a fourth order polyno- 
mial. Hence, the derivative dP/dV, and the bulk modulus 
B were calculated from the relation: B = —V(dP/dV), 
where V is the volume of the ZrC crystallographic unit 
cell. For LDA and GGA, the B values are given in Tabled 
and they agree quite well with experimental data. We also 
derived the bulk modulus from the energy-deformation 
relationship and the stress-strain relations, and the last 
proved to be more accurate than differentiation of the to- 
tal energy curve (see Tabled]). 



We used stress-strain relations to obtain the elastic 
constants. The calculations of en, C12, C44 elastic constants 
involved elongation and shear deformations. Deformations 
from 0.5% to 3% in length and from 1 to 5 degrees in an- 
gle are used. The results of small and large deformations 
are consistent. The deformed lattices have space groups 
Fm3m, I4/mmm and I/mmm, for bulk expansion, elon- 
gation along z, and shear modes, respectively. In all these 
deformed lattices all atoms remain in high-symmetry sites 
so that the relaxation of the internal degrees of freedom 
during the supercell deformations are not necessary. The 
calculated elastic constants are compared with the exper- 
imental data in Tabled Generally they are in good agree- 
ment with experiments. The C12, C44 elastic constants are 
better reproduced by the GGA, but also en fits better 
with GGA, when experimental data are extrapolated to 
the stoichiometric ZrC (third row in Table [2]) . 



In summary, we have shown that the ab initio calcu- 
lations of HF forces together with the direct method lead 
to a satisfactory description of phonons and elastic con- 
stants in the ZrC crystal. The phonon dispersion relations 
and the density of states fit well to neutron experimen- 
tal data. The GGA offered better phonon frequencies and 
elastic constants. 
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